Energy loss of atoms at metal surfaces due to electron-hole pair excitations: 
First-principles theory of 'chemicurrents' 
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A method is presented for calculating electron-hole pair excitation due to an incident atom or 
molecule interacting with a metal surface. Energy loss is described using an ah initio approach that 
obtains a position-dependent friction coefficient for an adsorbate moving near a metal surface from 
a total energy pseudopotential calculation. A semi-classical forced oscillator model is constructed, 
using the same friction coefficient description of the energy loss, to describe excitation of the electron 
gas due to the incident molecule. This approach is applied to H and D atoms incident on a Cu(lll) 
surface, and we obtain theoretical estimates of the 'chemicurrents' measured by Nienhaus et al 
[Phys. Rev. Lett. 82, 446 (1999)] for these atoms incident on the surface of a Schottky diode. 

PACS numbers: 73.20.Hb, 34.50.Dy, 68.43.-h 
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Considerable progress has been made in recent years 
in the theory of gas-surface interactions. This has been 
due to the paraUel developments of large-scale electronic 
structure calculations based on density functional theory, 
combined with multi-dimensional quantum and classical 
analysis of the dynamics [J, [3| • Despite these advances 
there remains one key area that is still largely unexplored 
and poorly understood; the process of energy dissipa- 
tion into substrate degrees of freedom. Although this is 
known to be of central importance in many situations , 
there exist few 'real' calculations to date for the energy 
loss to either phonons or electrons in the surface Q . 

In particular there have been a number of recent ex- 
periments that provide convincing evidence that energy 
dissipation by the creation of electron-hole (e-/i) pairs is 
a significant effect in gas-surface dynamics [j, Q ■ Of par- 
ticular interest for this Letter are the results of Nienhaus 
and co-workers 0, 0] who directly measured the hot elec- 
trons and holes created at Ag and Cu surfaces by the 
adsorption of thermal hydrogen and deuterium atoms in 
the form of a 'chemicurrent' in a Schottky diode. Here 
we report ab initio calculations of e-h pair creation for 
H/Cu(lll) and the resulting chemicurrents. 

The calculation proceeds in three stages. First, a stan- 
dard Kohn-Sham (KS) total energy calculation is carried 
out for a range of positions of the incident atom. Sec- 
ond, the resulting KS states and potentials are used to 
obtain a friction coefhcient for the motion of the atom via 
Time Dependent Density Functional Theory (TDDFT). 
Finally, a classical trajectory for the incident atom is cal- 
culated and a forced oscillator model (FOM) is used to 
obtain a semi-classical (classical atomic motion coupled 
to quantum metallic electrons) description of e-h pair 



creation. It is important to note that out theory refers 
to a nearly adiabatic process with multiple low energy 
excitations. This is in contrast to the truly nonadiabatic 
charge transfer models used to describe, for example, ex- 
oemission of electrons [1] . 

The initial total energy calculation provides the ground 
state properties of an interacting surface/atom system 
using a plane-wave basis, pseudopotentials and a super- 
cell geometry. This super-cell method has the advantage 
that it retains the continuous spectrum of one-electron 
excitations, unlike cluster models Q, which is important 
for the interactions considered here. The resulting KS 
eigenstates are used to evaluate the friction coefficient, 
r], associated with the motion of an atom at a chosen po- 
sition and in a direction of choice, h, using the method 
described in a previous paper lol | . T his applies TDDFT 

ll| , and in essence con- 



together with a quasi-static limit 
sists of the evaluation of the expression 

7] = 27rh ^ ( eFQ:|h.Vy|eFQ! ) 



(1) 



where a are supplementary quantum numbers for states 
on the Fermi surface, and h.VT^ is the rate of change 
of the self-consistent KS potential with atomic position. 
As described in [l3| , care must be taken to interpret the 
super-cell geometry correctly when evaluating Eq. ^ 
within a plane-wave basis, since the system of interest 
consists of the isolated motion of an atom at a surface, 
while the super-cell naturally describes a coherent lattice 
in motion. It is also important to perform a correct dis- 
cretisation of Eq. Il]) within the finite available sampling 
of k space ■ 

These first two steps provide the total energy and fric- 
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tion coefficient at positions of choice and thus it is pos- 
sible to obtain the classical trajectory of an atom near 
the surface. This provides the classical dynamics of the 
interaction such as the total energy loss and the critical 
initial kinetic energy, Cc, defined as the initial energy be- 
low which the atom is trapped in the surface well. For 
simplicity, only the case of a single atom incident perpen- 
dicular to the surface is considered here, and this trajec- 
tory is denoted by z{t). 

Electron-hole pair excitations can be further analysed 
by implementing a FOM ^] description of the energy loss 
to the electron gas. This provides a straightforward, if 
approximate, quantum calculation of the response of the 
electron gas to the time varying self-consistent potential 
resulting from the classical motion of the incident atom. 



Past work [ij, [13|, [ij] has shown that the probability 



density that an incident atom will excite an e-h pair of 
energy huj, given that a single e-h excitation event has 
occured, can be written as 



ao 



(eFa|h.Vy|eFa')i(i)e"'"*dt 



(2) 

where ao = Ps{uj)duj is the total number of e-h pairs 
created by multi-excitation events. Eq. ([2]) uses the 'Bo- 
son Approximation IF described by Schonhammer and 
Gunnarsson [13, [3l ■ This approximation is constructed 
by using adiabatic electronic states in place of the true 
states and assuming that the matrix elements vary suffi- 
ciently slowly with energy that they can be approximated 
by the values taken at the Fermi energy. All information 
about the system and its evolution is contained in Ps{lu), 
together with the assumption that individual e-h exci- 
tation events are uncorrelated. The probability that an 
incident atom will excite n e-h pairs is then given by 
e~°'"aQ/n\. This allows a multi-excitation expansion to 
be used to obtain the total energy loss of the incident 
atom, although this is not necessary for our purposes 

a. 

Although Eq. ([2]) may be evaluated by the same ap- 
proach used to obtain the friction coefficient, this requires 
the storage and interpolation of a large number of ma- 
trix elements. Due to the approximations already made 
to develop the theory to this point this effort does not 
seem justifiable. Instead, we assume that the dependence 
of the matrix elements in Eq. ([2]) on z and a, a' can be 
expressed in the separable form fa.a'g{z)- Physically, 
this implies that each e-h pair excitation experiences the 
same time-dependent force, but with different strengths. 
This leads to 
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ri3 {z)z{t)e 



*dt 



(3) 



where different spin terms are shown explicitly (subscript 
a). The same approach has been used by previous au- 
thors [3, Ei , lB| , and it can be shown that the mean 



excitation energy of the electron gas (at any instant) is 
the same for Ps{u!) in Eqs. ([2]) and ([3]). 

The FOM may been used to calculate various physi- 
cally measurable quantities, such as the sticking proba- 
bility and average energy loss of the incident atom. Here 
we are interested in the number of electrons excited above 
a given energy Cs = hujs (relative to the Fermi energy) as 
a result of the multiple e-h excitations. This quantity is 
given by 



and is of particular interest because it provides a direct 
estimate of the number of electrons made available for de- 
tection as a 'chemicurrent' in the experiments described 
by Nienhaus et al. 0]. 

We also note that Eq. ^ does not take into account 
any loss of coherence due to the finite lifetime of excited 
electrons. Preliminary estimates of this effect show that 
a finite lifetime is important for heavy (i.e. slow mov- 
ing) adsorbates, but for the systems considered here the 
chemicurrent changes by at most 40% [l7| . 

H/Cu(lll) and D/Cu(lll) Calculations of the inter- 
action energy, friction coefficient, dynamics and e-h pair 
creation have been carried out for both H and D atoms 
moving perpendicular to the surface above the top site 
of Cu(lll). This system has been chosen for simplic- 
ity (although it should be noted that spin polarisation is 
needed to obtain the correct behaviour for large atom- 
surface separations) and because of its relevance to the 
'chemicurrent' experiments of Nienhaus et al. Cal- 
culations have also been carried out for H/D above the 
hep hollow site, and for motion parallel to the surface - 
hep results are similar to those for the top site in terms 
of e-h pair creation, so here we report only the latter. 
Full results will be presented elsewhere (l7j . 

The surface is modelled by a five-layer slab, with a 
vacuum gap equivalent to another five layers. A 2 x 2 in- 
plane super-cell is used - tests show that the deformation 
potential (h.Vy) caused by the displacement of H atoms 
is well localised within this unit cell. A spin-polarised 
version of the PW91 functional is used for exchange- 
correlation effects (see 18|, and references therein), a 
TrouUier-Martins [19[ pseudopotential is used for Cu, and 
H is represented by a Coulomb potential. The plane- 
wave, pseudopotential code CASTEP is used to obtain 
the self-consistent potentials and KS states. Calculations 
are performed with a plane-wave cut-off of 830 eV, 54 k- 
points are included in the full surface Brillouin zone, and 
a Fermi surface broadening of 0.25 eV is used. Total 
energy calculations are performed for a range of heights 
between 1.0 and 4.0 A above the surface. 

Friction coefficients are evaluated for the same range 
of heights and for motion perpendicular to the surface, 
as described by Trail et al. [10|. The natural way to do 
calculations of this type is to assume spin-adiabaticity 
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FIG. 1: Potential energy (solid line) and friction coefficient 
for H atom in perpendicular motion above the top site of 
Cu(lll). Dashed line is friction resulting from constrained 
spin calculations, and dotted line for the spin degenerate case. 



and minimise the total energy of the system with respect 
to the magnetisation density, as well as the charge den- 
sity. Surprisingly, this leads to unphysical results. As 
the atom approaches the surface a phase transition from 
spin polarised to spin-degenerate occurs, as the H affinity 
level crosses the Fermi level 17, 20i, and h.VK-r becomes 



singular at this point. This results in a singularity in 
■q and it is straightforward to show that this divergence 
leads to infinite stopping power, i.e. incident particles of 
any energy are stopped at the singular height. Numerical 
calculations confirm the existence of the singularity at a 
height of 2.39 A above the surface \vi\ . 

We take our assumption of spin-adiabaticity to be the 
source of this unrealistic behaviour. An alternative ap- 
proach is to assume that the period of time that the atom 
spends near the critical height is not large enough for spin 
relaxation to occur, leading to a formulation where the 
net polarisation of the system is held fixed at its initial 
value. 

Figure[T]shows the resulting potential energy curve and 
friction coefficient for the constrained spin calculations. 
This friction differs significantly from the fully spin re- 
laxed friction for 2 < z < 3 A; a peak is still present 
as the affinity level crosses the Fermi energy, but it is 
no longer singular. The friction coefficient for a spin- 
degenerate calculation (expected to be inaccurate for z 
above the critical height) is also shown for comparison. 
In the remainder of this paper constrained spin data are 
used for all results unless explicitly stated otherwise, but 
a comparison of constrained with degenerate spin will al- 
low us to assess the significance of the friction in different 
regions of space to our final results. 

A standard Verlet integration is carried out to obtain 
the classical trajectories for an incident H or D atom. 



based on the data of Fig. [T] and for a range of initial 
kinetic energies. The classical trajectories show expected 
features, with a critical initial kinetic energy of — 166 
(115) meV for H (D), and an escaping atom spending 
~ 0.04 ps in the potential weU [z(t) < 3.0 A]. 

We now turn to the number of electrons excited above a 
given energy eg (relative to the Fermi energy), Ng[e > Es)- 
For a Schottky barrier of 0.6-0.65 eV, and a thermal flux 
of H atoms (300 - 350K) perpendicularly incident on the 
clean Cu surface of a Cu/Si Schottky diode at 135K, 
Nienhaus et al. Q detected 1.5 x 10~^ electrons/atom. 

Equation ([4]) is applied to the results of the FOM cal- 
culations, with = |/cT = 38.8 meV (T = 300K; the 
results vary slowly with incident energy in the thermal 
range, and taking a Boltzmann average makes no discern- 
able difference). Previous applications of the FOM have 
employed elastic trajectories, as this allows the calcula- 
tion of the probability that the incident atom will lose 
a given energy after traversing the adsorption well twice 
(i.e. in and out). The probability that this energy loss is 
greater than the incident energy of the atom gives a good 
estimate of the sticking probability, the quantity of inter- 
est in previous calculations. However, we are interested 
in all of the excitations that result from the atom relax- 
ing into the adsorption well, hence we use a full inelastic 
trajectory to obtain the chemicurrent. 

Results are shown for a range of Schottky barrier 
heights in Fig. [2l for H and D, and for both spin- 
constrained and spin-degenerate friction coefficients. In 
both cases the spin-constrained potential is used - the 
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FIG. 2: Number of electrons per atom made available for 
detection over a Schottky barrier of height es for H and D in- 
cident on the top site of Cu(lll). Solid line is for friction 
resulting from constrained spin, and dashed line for spin- 
degenerate. The vertical lines span the range of Schottky 
barrier heights found by Nienhaus et al. for both Cu and Ag. 
Inset shows the ratio of the chemicurrents due to H and D, 
for constrained spin. 
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spin-degenerate calculation results in an unrealistically 
deep potential. Spin-degenerate chemicurrents differ 
from the spin-constrained case by < 20% for the barrier 
heights of interest. The similarity of the results for very 
different friction functions is a consequence of the inci- 
dent atom spending a small proportion of the total tra- 
jectory in the region of space where the spin-degenerate 
and spin-constrained friction differ. If a large (though 
non-singular) friction is present around z = 2.39 A this 
will have little effect on the chemicurrent since the atom 
spends a short period of time and has a small velocity in 
this region (Eq. [S]). We conclude that a precise descrip- 
tion of friction in the region of the spin-phase transition 
is relatively unimportant, provided an unphysical singu- 
larity is avoided. 

Before any comparison with experimental results may 
be made, detection efficiency must be estimated. As dis- 
cussed by Nienhaus two effects dominate the attenua- 
tion of electron current in their Schottky detectors. First 
we consider the loss of ballistic electrons due to scatter- 
ing events in the thin metal film. This was found to 
be well described by an exponential decay, and for the 
device used to obtain the results above a ballistic mean 
free path (mfp) of 175 A was measured. For the film 
thickness of 75 A this results in attenuation by a factor 
of 0.65. Second, a transmission coefficient exists at the 
Cu/Si interface. This may be estimated from the results 
of Nienhaus 0] by comparison of the Cu/Si device with a 
multilayer device constructed by the deposition of a thin 
Cu film onto a Ag/Si detector, suggesting attenuation by 
a factor of 0.01 - 0.02. 

Correcting for these factors, and combining the un- 
certainties in the barrier heights and total attenuation 
factors, we may estimate a detected chemicurrent of 
(0.63 — 1.90) X 10~* electrons/atom, which is in good 
agreement with experiment. We may also obtain a rough 
estimate of the chemicurrent for a Ag Schottky detector 
by using the experimentally measured Ag barrier height 
(0.5-0.55 eV) in our Cu calculations. Correcting for an 
experimental mfp of 50 A [^0] (and Ag film thick- 
ness of 75 A) and using an estimated transmission of 
~ 1 at the Ag/Si interface results in a chemicurrent of 
(5.6 — 10.7) X 10"'^ electrons/atom for H. This again com- 
pares well with the experimental result of 4.5 x 10""^ elec- 
trons/atom |f|. 

An interesting signature of the experimental chemicur- 
rent is the strong isotope dependence j3l ■ Our calculation 
predicts an unusually strong isotope effect due to the es- 
sentially exponential variation of the chemicurrent with 
barrier height, and simple energy scaling in the dynam- 
ics (we note that scaling the energy variable by \/2 in 
iVe(e > Cs) results in exact agreement between the chemi- 
currents for H and D for an elastic trajectory, and this 
holds approximately for the inelastic case). Nienhaus et 
at reported no results for D incident on Cu, but for Ag 
the detected number of electrons/atom for H was larger 



than for D by a factor of about 6. The inset in Fig. [2] 
shows the calculated ratio of chemicurrents for H and D, 
and in the range of relevant barrier heights we find an 
isotope effect close to the experimental result. 

In summary, an ab initio description of the static ad- 
sorbate/metal surface and a dynamic FOM have been 
used to describe e-h pair creation. This approach in- 
cludes energy loss by the classical incident atom consis- 
tently, and our results for H/Cu(lll) are in good agree- 
ment with the chemicurrents measured by Nienhaus et al. 
Calculations show that simple dynamics can account for 
the large difference between the measured chemicurrents 
for H and D. 
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